home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / graphics / povray2.zip / SOURCE / TRIANGLE.C < prev    next >
C/C++ Source or Header  |  1993-10-25  |  20KB  |  695 lines
  1. /****************************************************************************
  2. *                triangle.c
  3. *
  4. *  This module implements primitives for triangles and smooth triangles.
  5. *
  6. *  from Persistence of Vision Raytracer
  7. *  Copyright 1993 Persistence of Vision Team
  8. *---------------------------------------------------------------------------
  9. *  NOTICE: This source code file is provided so that users may experiment
  10. *  with enhancements to POV-Ray and to port the software to platforms other 
  11. *  than those supported by the POV-Ray Team.  There are strict rules under
  12. *  which you are permitted to use this file.  The rules are in the file
  13. *  named POVLEGAL.DOC which should be distributed with this file. If 
  14. *  POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  15. *  Team Coordinator by leaving a message in CompuServe's Graphics Developer's
  16. *  Forum.  The latest version of POV-Ray may be found there as well.
  17. *
  18. * This program is based on the popular DKB raytracer version 2.12.
  19. * DKBTrace was originally written by David K. Buck.
  20. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  21. *
  22. *****************************************************************************/
  23.  
  24. #include "frame.h"
  25. #include "vector.h"
  26. #include "povproto.h"
  27.  
  28. METHODS Triangle_Methods = 
  29.   {
  30.   All_Triangle_Intersections,
  31.   Inside_Triangle, Triangle_Normal,
  32.   Copy_Triangle,
  33.   Translate_Triangle, Rotate_Triangle,
  34.   Scale_Triangle, Transform_Triangle, Invert_Triangle, Destroy_Triangle
  35.   };
  36.  
  37. METHODS Smooth_Triangle_Methods = 
  38.   {
  39.   All_Triangle_Intersections,
  40.   Inside_Triangle, Smooth_Triangle_Normal,
  41.   Copy_Smooth_Triangle,
  42.   Translate_Smooth_Triangle, Rotate_Smooth_Triangle,
  43.   Scale_Smooth_Triangle, Transform_Smooth_Triangle, 
  44.   Invert_Smooth_Triangle, Destroy_Triangle
  45.   };
  46.  
  47. extern RAY *CM_Ray;
  48. extern long Ray_Triangle_Tests, Ray_Triangle_Tests_Succeeded;
  49.  
  50. #define max3(x,y,z) ((x>y)?((x>z)?1:3):((y>z)?2:3))
  51.  
  52. #define MAX3(x,y,z) (((x)>(y))?(((x)>(z))?(x):(z)):(((y)>(z))?(y):(z)))
  53. #define MIN3(x,y,z) (((x)<(y))?(((x)<(z))?(x):(z)):(((y)<(z))?(y):(z)))
  54.  
  55. void Find_Triangle_Dominant_Axis(Triangle)
  56. TRIANGLE *Triangle;
  57.   {
  58.   DBL x, y, z;
  59.  
  60.   x = fabs(Triangle->Normal_Vector.x);
  61.   y = fabs (Triangle->Normal_Vector.y);
  62.   z = fabs (Triangle->Normal_Vector.z);
  63.   switch (max3(x, y, z)) 
  64.   {
  65.   case 1: 
  66.     Triangle->Dominant_Axis = X_AXIS;
  67.     break;
  68.   case 2: 
  69.     Triangle->Dominant_Axis = Y_AXIS;
  70.     break;
  71.   case 3: 
  72.     Triangle->Dominant_Axis = Z_AXIS;
  73.     break;
  74.   }
  75.   }
  76.  
  77. void Compute_Smooth_Triangle (Triangle)
  78. SMOOTH_TRIANGLE *Triangle;
  79.   {
  80.   VECTOR P3MinusP2, VTemp1, VTemp2;
  81.   DBL x, y, z, uDenominator, Proj;
  82.  
  83.   VSub (P3MinusP2, Triangle->P3, Triangle->P2);
  84.   x = fabs (P3MinusP2.x);
  85.   y = fabs (P3MinusP2.y);
  86.   z = fabs (P3MinusP2.z);
  87.  
  88.   switch (max3 (x, y, z)) 
  89.   {
  90.   case 1:  
  91.     Triangle->vAxis = X_AXIS;
  92.     Triangle->BaseDelta = P3MinusP2.x;
  93.     break;
  94.  
  95.   case 2:  
  96.     Triangle->vAxis = Y_AXIS;
  97.     Triangle->BaseDelta = P3MinusP2.y;
  98.     break;
  99.  
  100.   case 3:  
  101.     Triangle->vAxis = Z_AXIS;
  102.     Triangle->BaseDelta = P3MinusP2.z;
  103.     break;
  104.   }   
  105.  
  106.   VSub (VTemp1, Triangle->P2, Triangle->P3);
  107.   VNormalize (VTemp1, VTemp1);
  108.   VSub (VTemp2, Triangle->P1, Triangle->P3);
  109.   VDot (Proj, VTemp2, VTemp1);
  110.   VScaleEq (VTemp1, Proj);
  111.   VSub (Triangle->Perp, VTemp1, VTemp2);
  112.   VNormalize (Triangle->Perp, Triangle->Perp);
  113.   VDot (uDenominator, VTemp2, Triangle->Perp);
  114.   VInverseScaleEq (Triangle->Perp, -uDenominator);
  115.   }
  116.  
  117. int Compute_Triangle (Triangle,Smooth)
  118. TRIANGLE *Triangle;
  119. int Smooth;
  120.   {
  121.   VECTOR V1, V2, Temp;
  122.   DBL Length, T1, T2, T3;
  123.  
  124.   VSub (V1, Triangle->P1, Triangle->P2);
  125.   VSub (V2, Triangle->P3, Triangle->P2);
  126.   VCross (Triangle->Normal_Vector, V1, V2);
  127.   VLength (Length, Triangle->Normal_Vector);
  128.   /* Set up a flag so we can ignore degenerate triangles */
  129.   if (Length < 1.0e-9)
  130.     {
  131.     Triangle->Degenerate_Flag = TRUE;
  132.     return (0);
  133.     }
  134.  
  135.   /* Normalize the normal vector. */
  136.   VInverseScaleEq (Triangle->Normal_Vector, Length);
  137.  
  138.   VDot (Triangle->Distance, Triangle->Normal_Vector, Triangle->P1);
  139.   Triangle->Distance *= -1.0;
  140.   Find_Triangle_Dominant_Axis(Triangle);
  141.  
  142.   switch (Triangle->Dominant_Axis) 
  143.   {
  144.   case X_AXIS:
  145.     if ((Triangle->P2.y - Triangle->P3.y)*(Triangle->P2.z - Triangle->P1.z) <
  146.       (Triangle->P2.z - Triangle->P3.z)*(Triangle->P2.y - Triangle->P1.y)) 
  147.       {
  148.       Temp = Triangle->P2;
  149.       Triangle->P2 = Triangle->P1;
  150.       Triangle->P1 = Temp;
  151.       if (Smooth) 
  152.         {
  153.         Temp = ((SMOOTH_TRIANGLE *) Triangle)->N2;
  154.         ((SMOOTH_TRIANGLE *) Triangle)->N2 = ((SMOOTH_TRIANGLE *) Triangle)->N1;
  155.         ((SMOOTH_TRIANGLE *) Triangle)->N1 = Temp;
  156.         }
  157.       }
  158.     break;
  159.  
  160.   case Y_AXIS:
  161.     if ((Triangle->P2.x - Triangle->P3.x)*(Triangle->P2.z - Triangle->P1.z) <
  162.       (Triangle->P2.z - Triangle->P3.z)*(Triangle->P2.x - Triangle->P1.x)) 
  163.       {
  164.       Temp = Triangle->P2;
  165.       Triangle->P2 = Triangle->P1;
  166.       Triangle->P1 = Temp;
  167.       if (Smooth) 
  168.         {
  169.         Temp = ((SMOOTH_TRIANGLE *) Triangle)->N2;
  170.         ((SMOOTH_TRIANGLE *) Triangle)->N2 = ((SMOOTH_TRIANGLE *) Triangle)->N1;
  171.         ((SMOOTH_TRIANGLE *) Triangle)->N1 = Temp;
  172.         }
  173.       }
  174.     break;
  175.  
  176.   case Z_AXIS:
  177.     if ((Triangle->P2.x - Triangle->P3.x)*(Triangle->P2.y - Triangle->P1.y) <
  178.       (Triangle->P2.y - Triangle->P3.y)*(Triangle->P2.x - Triangle->P1.x)) 
  179.       {
  180.       Temp = Triangle->P2;
  181.       Triangle->P2 = Triangle->P1;
  182.       Triangle->P1 = Temp;
  183.       if (Smooth) 
  184.         {
  185.         Temp = ((SMOOTH_TRIANGLE *) Triangle)->N2;
  186.         ((SMOOTH_TRIANGLE *) Triangle)->N2 = ((SMOOTH_TRIANGLE *) Triangle)->N1;
  187.         ((SMOOTH_TRIANGLE *) Triangle)->N1 = Temp;
  188.         }
  189.       }
  190.     break;
  191.   }
  192.  
  193.   if (Smooth)
  194.     Compute_Smooth_Triangle((SMOOTH_TRIANGLE *) Triangle);
  195.  
  196.   /* Build the bounding information from the vertices */
  197.   /* Use temps so macro not too big. */
  198.   T1 = MIN3(Triangle->P1.x, Triangle->P2.x, Triangle->P3.x);
  199.   T2 = MIN3(Triangle->P1.y, Triangle->P2.y, Triangle->P3.y);
  200.   T3 = MIN3(Triangle->P1.z, Triangle->P2.z, Triangle->P3.z);
  201.   Make_Vector(&Triangle->Bounds.Lower_Left,T1,T2,T3);
  202.  
  203.   T1 = MAX3(Triangle->P1.x, Triangle->P2.x, Triangle->P3.x);
  204.   T2 = MAX3(Triangle->P1.y, Triangle->P2.y, Triangle->P3.y);
  205.   T3 = MAX3(Triangle->P1.z, Triangle->P2.z, Triangle->P3.z);
  206.   Make_Vector(&Triangle->Bounds.Lengths,T1,T2,T3);
  207.  
  208.   VSub(Triangle->Bounds.Lengths, Triangle->Bounds.Lengths,
  209.     Triangle->Bounds.Lower_Left);
  210.   Triangle->Bounds.Lengths.x += EPSILON;
  211.   Triangle->Bounds.Lengths.y += EPSILON;
  212.   Triangle->Bounds.Lengths.z += EPSILON;
  213.   return (1);
  214.   }
  215.  
  216. int All_Triangle_Intersections (Object, Ray, Depth_Stack)
  217. OBJECT *Object;
  218. RAY *Ray;
  219. ISTACK *Depth_Stack;
  220.   {
  221.   DBL Depth;
  222.   VECTOR IPoint;
  223.  
  224.   if (Intersect_Triangle (Ray, (TRIANGLE *)Object, &Depth))
  225.     {
  226.     VScale (IPoint, Ray -> Direction, Depth);
  227.     VAddEq (IPoint, Ray -> Initial);
  228.     if (Point_In_Clip(&IPoint,Object->Clip))
  229.       {
  230.       push_entry(Depth,IPoint,Object,Depth_Stack);
  231.       return (TRUE);
  232.       }
  233.     }
  234.   return (FALSE);
  235.   }
  236.  
  237. int Intersect_Triangle (Ray, Triangle, Depth)
  238. RAY *Ray;
  239. TRIANGLE *Triangle;
  240. DBL *Depth;
  241.   {
  242.   DBL NormalDotOrigin, NormalDotDirection;
  243.   DBL s, t;
  244.  
  245.   Ray_Triangle_Tests++;
  246.   if(Triangle->Degenerate_Flag)
  247.     return(FALSE);
  248.  
  249.   if (Ray == CM_Ray) 
  250.     {
  251.     if (!Triangle->CMCached) 
  252.       {
  253.       VDot (Triangle->CMNormDotOrigin, Triangle->Normal_Vector, Ray->Initial);
  254.       Triangle->CMNormDotOrigin += Triangle->Distance;
  255.       Triangle->CMNormDotOrigin *= -1.0;
  256.       Triangle->CMCached = TRUE;
  257.       }
  258.  
  259.     VDot (NormalDotDirection, Triangle->Normal_Vector, Ray->Direction);
  260.     if ((NormalDotDirection < Small_Tolerance) &&
  261.       (NormalDotDirection > -Small_Tolerance))
  262.       return (FALSE);
  263.  
  264.     *Depth = Triangle->CMNormDotOrigin / NormalDotDirection;
  265.     }
  266.   else 
  267.     {
  268.     VDot (NormalDotOrigin, Triangle->Normal_Vector, Ray->Initial);
  269.     NormalDotOrigin += Triangle->Distance;
  270.     NormalDotOrigin *= -1.0;
  271.  
  272.     VDot (NormalDotDirection, Triangle->Normal_Vector, Ray->Direction);
  273.     if ((NormalDotDirection < Small_Tolerance) &&
  274.       (NormalDotDirection > -Small_Tolerance))
  275.       return (FALSE);
  276.  
  277.     *Depth = NormalDotOrigin / NormalDotDirection;
  278.     }
  279.  
  280.   if ((*Depth < Small_Tolerance) || (*Depth > Max_Distance))
  281.     return (FALSE);
  282.  
  283.   switch (Triangle->Dominant_Axis) 
  284.   {
  285.   case X_AXIS:
  286.     s = Ray->Initial.y + *Depth * Ray->Direction.y;
  287.     t = Ray->Initial.z + *Depth * Ray->Direction.z;
  288.  
  289.     if ((Triangle->P2.y - s)*(Triangle->P2.z - Triangle->P1.z) <
  290.       (Triangle->P2.z - t)*(Triangle->P2.y - Triangle->P1.y))
  291.       return (FALSE);
  292.  
  293.     if ((Triangle->P3.y - s)*(Triangle->P3.z - Triangle->P2.z) <
  294.       (Triangle->P3.z - t)*(Triangle->P3.y - Triangle->P2.y))
  295.       return (FALSE);
  296.  
  297.     if ((Triangle->P1.y - s)*(Triangle->P1.z - Triangle->P3.z) <
  298.       (Triangle->P1.z - t)*(Triangle->P1.y - Triangle->P3.y))
  299.       return (FALSE);
  300.  
  301.     Ray_Triangle_Tests_Succeeded++;
  302.     return (TRUE);
  303.  
  304.   case Y_AXIS:
  305.     s = Ray->Initial.x + *Depth * Ray->Direction.x;
  306.     t = Ray->Initial.z + *Depth * Ray->Direction.z;
  307.  
  308.     if ((Triangle->P2.x - s)*(Triangle->P2.z - Triangle->P1.z) <
  309.       (Triangle->P2.z - t)*(Triangle->P2.x - Triangle->P1.x))
  310.       return (FALSE);
  311.  
  312.     if ((Triangle->P3.x - s)*(Triangle->P3.z - Triangle->P2.z) <
  313.       (Triangle->P3.z - t)*(Triangle->P3.x - Triangle->P2.x))
  314.       return (FALSE);
  315.  
  316.     if ((Triangle->P1.x - s)*(Triangle->P1.z - Triangle->P3.z) <
  317.       (Triangle->P1.z - t)*(Triangle->P1.x - Triangle->P3.x))
  318.       return (FALSE);
  319.  
  320.     Ray_Triangle_Tests_Succeeded++;
  321.     return (TRUE);
  322.  
  323.   case Z_AXIS:
  324.     s = Ray->Initial.x + *Depth * Ray->Direction.x;
  325.     t = Ray->Initial.y + *Depth * Ray->Direction.y;
  326.  
  327.     if ((Triangle->P2.x - s)*(Triangle->P2.y - Triangle->P1.y) <
  328.       (Triangle->P2.y - t)*(Triangle->P2.x - Triangle->P1.x))
  329.       return (FALSE);
  330.  
  331.     if ((Triangle->P3.x - s)*(Triangle->P3.y - Triangle->P2.y) <
  332.       (Triangle->P3.y - t)*(Triangle->P3.x - Triangle->P2.x))
  333.       return (FALSE);
  334.  
  335.     if ((Triangle->P1.x - s)*(Triangle->P1.y - Triangle->P3.y) <
  336.       (Triangle->P1.y - t)*(Triangle->P1.x - Triangle->P3.x))
  337.       return (FALSE);
  338.  
  339.     Ray_Triangle_Tests_Succeeded++;
  340.     return (TRUE);
  341.   }
  342.   return (FALSE);
  343.   }
  344.  
  345. int Inside_Triangle (IPoint, Object)
  346. VECTOR *IPoint;
  347. OBJECT *Object;
  348.   {
  349.   return (FALSE);
  350.   }
  351.  
  352. void Triangle_Normal (Result, Object, IPoint)
  353. OBJECT *Object;
  354. VECTOR *Result, *IPoint;
  355.   {
  356.   *Result = ((TRIANGLE *)Object)->Normal_Vector;
  357.   }
  358.  
  359. void *Copy_Triangle (Object)
  360. OBJECT *Object;
  361.   {
  362.   TRIANGLE *New;
  363.  
  364.   New = Create_Triangle ();
  365.   *New = * ((TRIANGLE *)Object);
  366.  
  367.   return (New);
  368.   }
  369.  
  370. void Translate_Triangle (Object, Vector)
  371. OBJECT *Object;
  372. VECTOR *Vector;
  373.   {
  374.   TRIANGLE *Triangle = (TRIANGLE *) Object;
  375.   VECTOR Translation;
  376.  
  377.   if(Triangle->Degenerate_Flag) return;
  378.  
  379.   VEvaluate (Translation, Triangle->Normal_Vector, *Vector);
  380.   Triangle->Distance -= Translation.x + Translation.y + Translation.z;
  381.   VAddEq (Triangle->P1, *Vector)
  382.     VAddEq (Triangle->P2, *Vector)
  383.       VAddEq (Triangle->P3, *Vector)
  384.  
  385.         /* Recalculate the bounds */
  386.         VAddEq(Object->Bounds.Lower_Left, *Vector);
  387.   }
  388.  
  389. void Rotate_Triangle (Object, Vector)
  390. OBJECT *Object;
  391. VECTOR *Vector;
  392.   {
  393.   TRANSFORM Trans;
  394.  
  395.   if(((TRIANGLE *)Object)->Degenerate_Flag) return;
  396.  
  397.   Compute_Rotation_Transform (&Trans, Vector);
  398.   Transform_Triangle (Object, &Trans);
  399.   }
  400.  
  401. void Scale_Triangle (Object, Vector)
  402. OBJECT *Object;
  403. VECTOR *Vector;
  404.   {
  405.   TRIANGLE *Triangle = (TRIANGLE *) Object;
  406.   DBL Length,T1,T2,T3;
  407.  
  408.   if(Triangle->Degenerate_Flag) return;
  409.  
  410.   Triangle->Normal_Vector.x = Triangle->Normal_Vector.x / Vector->x;
  411.   Triangle->Normal_Vector.y = Triangle->Normal_Vector.y / Vector->y;
  412.   Triangle->Normal_Vector.z = Triangle->Normal_Vector.z / Vector->z;
  413.  
  414.   VLength(Length, Triangle->Normal_Vector);
  415.   VInverseScaleEq (Triangle->Normal_Vector, Length);
  416.   Triangle->Distance /= Length;
  417.  
  418.   VEvaluateEq (Triangle->P1, *Vector);
  419.   VEvaluateEq (Triangle->P2, *Vector);
  420.   VEvaluateEq (Triangle->P3, *Vector);
  421.  
  422.   /* Recompute the bounds */
  423.   /* Use temps so macro not too big. */
  424.   T1 = MIN3(Triangle->P1.x, Triangle->P2.x, Triangle->P3.x);
  425.   T2 = MIN3(Triangle->P1.y, Triangle->P2.y, Triangle->P3.y);
  426.   T3 = MIN3(Triangle->P1.z, Triangle->P2.z, Triangle->P3.z);
  427.   Make_Vector(&Triangle->Bounds.Lower_Left,T1,T2,T3);
  428.  
  429.   T1 = MAX3(Triangle->P1.x, Triangle->P2.x, Triangle->P3.x);
  430.   T2 = MAX3(Triangle->P1.y, Triangle->P2.y, Triangle->P3.y);
  431.   T3 = MAX3(Triangle->P1.z, Triangle->P2.z, Triangle->P3.z);
  432.   Make_Vector(&Triangle->Bounds.Lengths,T1,T2,T3);
  433.  
  434.   VSub(Triangle->Bounds.Lengths, Triangle->Bounds.Lengths,
  435.     Triangle->Bounds.Lower_Left);
  436.   Triangle->Bounds.Lengths.x += EPSILON;
  437.   Triangle->Bounds.Lengths.y += EPSILON;
  438.   Triangle->Bounds.Lengths.z += EPSILON;
  439.   }
  440.  
  441. void Transform_Triangle (Object, Trans)
  442. OBJECT *Object;
  443. TRANSFORM *Trans;
  444.   {
  445.   TRIANGLE *Triangle = (TRIANGLE *) Object;
  446.  
  447.   if(Triangle->Degenerate_Flag) return;
  448.  
  449.   MTransPoint (&Triangle->Normal_Vector,
  450.     &Triangle->Normal_Vector, Trans);
  451.   MTransPoint (&Triangle->P1, &Triangle->P1, Trans);
  452.   MTransPoint (&Triangle->P2, &Triangle->P2, Trans);
  453.   MTransPoint (&Triangle->P3, &Triangle->P3, Trans);
  454.   Compute_Triangle (Triangle,FALSE);
  455.   }
  456.  
  457. TRIANGLE *Create_Triangle()
  458.   {
  459.   TRIANGLE *New;
  460.  
  461.   if ((New = (TRIANGLE *) malloc (sizeof (TRIANGLE))) == NULL)
  462.     MAError ("triangle");
  463.  
  464.   INIT_OBJECT_FIELDS(New,TRIANGLE_OBJECT,&Triangle_Methods)
  465.  
  466.     Make_Vector (&(New -> Normal_Vector), 0.0, 1.0, 0.0);
  467.   New -> Distance = 0.0;
  468.   New -> CMNormDotOrigin = 0.0;
  469.   New -> CMCached = FALSE;
  470.   Make_Vector (&(New -> P1), 0.0, 0.0, 0.0);
  471.   Make_Vector (&(New -> P2), 1.0, 0.0, 0.0);
  472.   Make_Vector (&(New -> P3), 0.0, 1.0, 0.0);
  473.   New -> Degenerate_Flag = FALSE;
  474.  
  475.   /* NOTE: Dominant_Axis is computed when Parse_Triangle calls
  476.    Compute_Triangle.  vAxis is used only for smooth triangles */
  477.  
  478.   return (New);
  479.   }
  480.  
  481. void Invert_Triangle (Object)
  482. OBJECT *Object;
  483.   {
  484.   return;
  485.   }
  486.  
  487. /* Calculate the Phong-interpolated vector within the triangle
  488.    at the given intersection point. The math for this is a bit
  489.    bizarre:
  490.  
  491.     -         P1
  492.     |        /|\ \
  493.     |       / |Perp\
  494.     |      /  V  \   \
  495.     |     /   |    \   \
  496.   u |    /____|_____PI___\
  497.     |   /     |       \    \
  498.     -  P2-----|--------|----P3
  499.               Pbase    PIntersect
  500.         |-------------------|
  501.                        v
  502.  
  503.    Triangle->Perp is a unit vector from P1 to Pbase. We calculate
  504.  
  505.    u = (PI - P1) DOT Perp / ((P3 - P1) DOT Perp).
  506.  
  507.    We then calculate where the line from P1 to PI intersects the line P2 to P3:
  508.    PIntersect = (PI - P1)/u.
  509.  
  510.    We really only need one coordinate of PIntersect.  We then calculate v as:
  511.  
  512.       v = PIntersect.x / (P3.x - P2.x)
  513.  or   v = PIntersect.y / (P3.y - P2.y)
  514.  or   v = PIntersect.z / (P3.z - P2.z)
  515.  
  516.    depending on which calculation will give us the best answers.
  517.  
  518.    Once we have u and v, we can perform the normal interpolation as:
  519.  
  520.      NTemp1 = N1 + u(N2 - N1);
  521.      NTemp2 = N1 + u(N3 - N1);
  522.      Result = normalize (NTemp1 + v(NTemp2 - NTemp1))
  523.  
  524.    As always, any values which are constant for the triangle are cached
  525.    in the triangle.
  526. */
  527.  
  528. void Smooth_Triangle_Normal (Result, Object, IPoint)
  529. OBJECT *Object;
  530. VECTOR *Result, *IPoint;
  531.   {
  532.   SMOOTH_TRIANGLE *Triangle = (SMOOTH_TRIANGLE *) Object;
  533.   VECTOR PIMinusP1, NTemp1, NTemp2;
  534.   DBL u = 0.0, v = 0.0;
  535.  
  536.   VSub (PIMinusP1, *IPoint, Triangle->P1);
  537.   VDot (u, PIMinusP1, Triangle->Perp);
  538.   if (u < 1.0e-9) 
  539.     {
  540.     *Result = Triangle->N1;
  541.     return;
  542.     }
  543.  
  544.   /* BaseDelta contains P3.x-P2.x,  P3.y-P2.y, or P3.z-P2.z depending on the
  545.       value of vAxis. */
  546.  
  547.   switch (Triangle->vAxis) 
  548.   {
  549.   case X_AXIS:  
  550.     v = (PIMinusP1.x/u + Triangle->P1.x - Triangle->P2.x) / Triangle->BaseDelta;
  551.     break;
  552.  
  553.   case Y_AXIS:  
  554.     v = (PIMinusP1.y/u + Triangle->P1.y - Triangle->P2.y) / Triangle->BaseDelta;
  555.     break;
  556.  
  557.   case Z_AXIS:  
  558.     v = (PIMinusP1.z/u + Triangle->P1.z - Triangle->P2.z)/ Triangle->BaseDelta;
  559.     break;
  560.   }
  561.  
  562.   VSub (NTemp1, Triangle->N2, Triangle->N1);
  563.   VScaleEq (NTemp1, u);
  564.   VAddEq (NTemp1, Triangle->N1);
  565.   VSub (NTemp2, Triangle->N3, Triangle->N1);
  566.   VScaleEq (NTemp2, u);
  567.   VAddEq (NTemp2, Triangle->N1);
  568.   VSub (*Result, NTemp2, NTemp1);
  569.   VScaleEq (*Result, v);
  570.   VAddEq (*Result, NTemp1); 
  571.   VNormalize (*Result, *Result);
  572.   }
  573.  
  574. void *Copy_Smooth_Triangle (Object)
  575. OBJECT *Object;
  576.   {
  577.   SMOOTH_TRIANGLE *New;
  578.  
  579.   New = Create_Smooth_Triangle ();
  580.   *New = * ((SMOOTH_TRIANGLE *)Object);
  581.  
  582.   return (New);
  583.   }
  584.  
  585. void Translate_Smooth_Triangle (Object, Vector)
  586. OBJECT *Object;
  587. VECTOR *Vector;
  588.   {
  589.   SMOOTH_TRIANGLE *Triangle = (SMOOTH_TRIANGLE *) Object;
  590.   VECTOR Translation;
  591.  
  592.   if(Triangle->Degenerate_Flag) return;
  593.  
  594.   VEvaluate (Translation, Triangle->Normal_Vector, *Vector);
  595.   Triangle->Distance -= Translation.x + Translation.y + Translation.z;
  596.   VAddEq (Triangle->P1, *Vector)
  597.     VAddEq (Triangle->P2, *Vector)
  598.       VAddEq (Triangle->P3, *Vector)
  599.         Compute_Triangle ((TRIANGLE *) Triangle,TRUE);
  600.   }
  601.  
  602. void Rotate_Smooth_Triangle (Object, Vector)
  603. OBJECT *Object;
  604. VECTOR *Vector;
  605.   {
  606.   TRANSFORM Trans;
  607.  
  608.   if(((SMOOTH_TRIANGLE *)Object)->Degenerate_Flag) return;
  609.  
  610.   Compute_Rotation_Transform (&Trans, Vector);
  611.   Transform_Smooth_Triangle (Object, &Trans);
  612.   }
  613.  
  614. void Scale_Smooth_Triangle (Object, Vector)
  615. OBJECT *Object;
  616. VECTOR *Vector;
  617.   {
  618.   SMOOTH_TRIANGLE *Triangle = (SMOOTH_TRIANGLE *) Object;
  619.   DBL Length;
  620.  
  621.   if(Triangle->Degenerate_Flag) return;
  622.  
  623.   Triangle->Normal_Vector.x = Triangle->Normal_Vector.x / Vector->x;
  624.   Triangle->Normal_Vector.y = Triangle->Normal_Vector.y / Vector->y;
  625.   Triangle->Normal_Vector.z = Triangle->Normal_Vector.z / Vector->z;
  626.  
  627.   VLength(Length, Triangle->Normal_Vector);
  628.   VScaleEq (Triangle->Normal_Vector, 1.0 / Length);
  629.   Triangle->Distance /= Length;
  630.  
  631.   VEvaluateEq (Triangle->P1, *Vector);
  632.   VEvaluateEq (Triangle->P2, *Vector);
  633.   VEvaluateEq (Triangle->P3, *Vector);
  634.   Compute_Triangle ((TRIANGLE *) Triangle,TRUE);
  635.   }
  636.  
  637. void Transform_Smooth_Triangle (Object, Trans)
  638. OBJECT *Object;
  639. TRANSFORM *Trans;
  640.   {
  641.   SMOOTH_TRIANGLE *Triangle = (SMOOTH_TRIANGLE *) Object;
  642.  
  643.   if(Triangle->Degenerate_Flag) return;
  644.  
  645.   MTransPoint (&Triangle->Normal_Vector,
  646.     &Triangle->Normal_Vector, Trans);
  647.   MTransPoint (&Triangle->P1, &Triangle->P1, Trans);
  648.   MTransPoint (&Triangle->P2, &Triangle->P2, Trans);
  649.   MTransPoint (&Triangle->P3, &Triangle->P3, Trans);
  650.   MTransPoint (&Triangle->N1, &Triangle->N1, Trans);
  651.   MTransPoint (&Triangle->N2, &Triangle->N2, Trans);
  652.   MTransPoint (&Triangle->N3, &Triangle->N3, Trans);
  653.   Compute_Triangle ((TRIANGLE *) Triangle,TRUE);
  654.   }
  655.  
  656. void Invert_Smooth_Triangle (Object)
  657. OBJECT *Object;
  658.   {
  659.   return;
  660.   }
  661.  
  662. SMOOTH_TRIANGLE *Create_Smooth_Triangle()
  663.   {
  664.   SMOOTH_TRIANGLE *New;
  665.  
  666.   if ((New = (SMOOTH_TRIANGLE *) malloc (sizeof (SMOOTH_TRIANGLE))) == NULL)
  667.     MAError ("smooth triangle");
  668.  
  669.   INIT_OBJECT_FIELDS(New,SMOOTH_TRIANGLE_OBJECT,&Smooth_Triangle_Methods)
  670.  
  671.     Make_Vector (&(New->Normal_Vector), 0.0, 1.0, 0.0);
  672.   New->Distance = 0.0;
  673.   New -> CMNormDotOrigin = 0.0;
  674.   New -> CMCached = FALSE;
  675.   Make_Vector (&(New -> P1), 0.0, 0.0, 0.0);
  676.   Make_Vector (&(New -> P2), 1.0, 0.0, 0.0);
  677.   Make_Vector (&(New -> P3), 0.0, 1.0, 0.0);
  678.   Make_Vector (&(New -> N1), 0.0, 1.0, 0.0);
  679.   Make_Vector (&(New -> N2), 0.0, 1.0, 0.0);
  680.   Make_Vector (&(New -> N3), 0.0, 1.0, 0.0);
  681.   New -> BaseDelta = 0.0;
  682.   New -> Degenerate_Flag = FALSE;
  683.  
  684.   /* NOTE: Dominant_Axis and vAxis are computed when Parse_Triangle calls
  685.    Compute_Triangle.  */
  686.  
  687.   return (New);
  688.   }
  689.  
  690. void Destroy_Triangle (Object)
  691. OBJECT *Object;
  692.   {
  693.   free (Object);
  694.   }
  695.